#include <fstream>
#include <iostream>
#include <iomanip>
using namespace std;


#include "plant_1DSOS_RDM.h"
/*****************************************************************************/
plant_1DSOS_RDM::plant_1DSOS_RDM(double idt, int iL, int iMode, double iE_s, double iE_N){
	// Options
	LatticeSize = iL;
	mode = iMode;
	E_D = iE_D;
	E_S = iE_S;
	E_N = iE_N;

	// Inputs & outputs
	int input_num = 2;
	int input_dim[2] = {1,1};
	int output_num = 1;
	int output_dim[4] = {1,1,1,1};
	// Output #0: h
	// Output #1: r2
	// Output #2: a
	// Output #3: m2
	block::initialize(input_num,input_dim,output_num,output_dim,idt);
	
	// Monte Carlo Simulators
	MC_Event = new MCSimulator(2);    // 0 = adsorption, 1 = migration, 2 = desorption
	MC_AdsorptionSite = new MCSimulator(LatticeSize);
	MC_MigrationSite = new MCSimulator(LatticeSize);
	MC_Neighbor = new MCSimulator(2); // 0 = left, 1 = right

	int i;
	// All site have the same probability of adsorption
	for(i=0;i<LatticeSize;i++){
		MC_AdsorptionSite->UpdateProb(i,1);
	}
	// The probability list for migration is constructed in reset()

	// Right and left have the same probability of migration to
	for(i=0;i<2;i++){
		MC_Neighbor->UpdateProb(i,1);
	}

}
/*****************************************************************************/
plant_1DSOS_RDM::~plant_1DSOS_RDM(){
	delete MC_Event;
	delete MC_AdsorptionSite;
	delete MC_MigrationSite;
	delete MC_Neighbor;
}
/*****************************************************************************/
void plant_1DSOS_RDM::reset(){
	
}
/*****************************************************************************/
void plant_1DSOS_RDM::update(double sysTime){
	int EventID = MC_Event->TryOnce();
	switch EventID{
		case 0:
			{adsorption();break;}
		case 1:
			{migration();break;}
		case 2:
			{desorption();break;}
	}
}
/*****************************************************************************/
void plant_1DSOS_RDM::adsorption(){
	int SiteID = MC_AdsorptionSite->TryOnce();
	add_remove_atom(SiteID,1);
}
/*****************************************************************************/
void plant_1DSOS_RDM::migration(){
	int SiteID = MC_MigrationSite->TryOnce();
	int NeighborID = MC_Neighbor->TryOnce();
	add_remove_atom(SiteID,-1);
}
/*****************************************************************************/
void plant_1DSOS_RDM::init_migr_rate()
{
	int i;
	double k_m;

	k_m=2*temperature*KB2H;

	for(i=0;i<N_CL-1;i++)
	{
		MIGR_RT[i]=k_m*exp(-(E_S+i*E_N)*EV2KB/(temperature+0.0));
	}
    MIGR_RT[N_CL-1]=0;
	return;
}
/*****************************************************************************/
void plant_1DSOS_RDM::init_desorp_rate()
{
	int i;
	double k_m;

	k_m=2*temperature*KB2H;

	for(i=0;i<N_CL;i++)
	{
		//DESORP_RT[i]=k_m*exp(-(E_D+i*E_N)*EV2KB/(temperature+0.0));
		DESORP_RT[i]=0;
	}
	return;
}